package edu.uaskl.jqf.tools.matrix.solvers;

import java.util.ArrayList;
import java.util.List;

import edu.uaskl.jqf.model.Complex;
import edu.uaskl.jqf.model.matrix.MatrixComplex;
import edu.uaskl.jqf.model.matrix.MatrixComplexQuadraticIdentity;
import edu.uaskl.jqf.model.matrix.MatrixSolverEigen;
import edu.uaskl.jqf.model.matrix.MatrixTransformationHessenberg;
import edu.uaskl.jqf.tools.ArrayTools;

public enum Eigen {
    /** Singleton */
    INSTANCE;

    private static final boolean IS_PRINT_ACTIVATED = false;
    private static final double ACCURACY = Math.pow(2.0, -23);

    public MatrixSolverEigen solve(final MatrixTransformationHessenberg matrixTransformationHessenberg) {
        if (matrixTransformationHessenberg.getMatrixComplex().getRows() == 2) // then it must be 2x2, because hessenberg is quadratic
            return getEigenValuesFrom2x2MatrixComplex(matrixTransformationHessenberg.getMatrixComplex());

        // based on http://www.csee.umbc.edu/~squire/cs455_l13.html by Jon Squire
        final Complex[][] A = matrixTransformationHessenberg.getMatrixComplex().getValuesCopy();
        final int[] rowcol = matrixTransformationHessenberg.getTransformationVectorCopy();

        double anorm = 0.0;
        final int size = A.length;

        print("acc=" + ACCURACY + " = 2^-23");
        Complex T = Complex.get00();
        int sizeMul30 = size * 30; // heuristic on maximum iterations
        final Complex[] eigenValues = new Complex[size];
        final Complex[][] eigenVectors = new MatrixComplexQuadraticIdentity(size).getValuesCopy();

        // starting from Hessenberg reduction
        // for i in reverse A'FIRST+1..A'LAST-1 loop
        for (int ii = size - 2; ii > 0; ii--) {
            final int j = rowcol[ii];
            for (int k = ii + 1; k < size; k++)
                eigenVectors[k][ii] = A[k][ii - 1];
            if (ii != j) {
                // for k in i..A'LAST loop
                for (int k = ii; k < size; k++) {
                    eigenVectors[ii][k] = eigenVectors[j][k];
                    eigenVectors[j][k] = Complex.get00();
                }
                eigenVectors[j][ii] = Complex.get10();
            }
        }
        int len = size - 1; // used as subscript, loop test <=ien

        // ien is decremented
        while (len >= 0) { // 260
            print("260 low=0, ien=" + len);
            int its = 0;
            // look for small single subdiagonal element
            while (true) {
                print("in 280");
                int k = 0;
                // for kk in reverse low+1..ien loop // 300
                for (int kk = len; kk > 0; kk--) {
                    print("300 kk=" + kk);

                    if (A[kk][kk - 1].getSumAbs() <= (ACCURACY * (A[kk - 1][kk - 1].getSumAbs() + A[kk][kk].getSumAbs()))) {
                        k = kk;
                        break;
                    }
                } // 300
                print("exiting 300 with k=" + k);
                if (k == len)
                    break; // exit while true loop
                if (sizeMul30 <= 0)
                    throw new IllegalStateException("Maximum iteration count reached, exit. count: " + (size * 30));

                // compute shift
                Complex S = null;
                if ((its == 10) || (its == 20))
                    S = new Complex(Math.abs(A[len][len - 1].getRealPart()) + Math.abs(A[len - 1][len - 2].getRealPart()), Math.abs(A[len][len - 1]
                            .getImaginaryPart()) + Math.abs(A[len - 1][len - 2].getImaginaryPart()));
                else {
                    S = A[len][len];
                    final Complex x = A[len - 1][len].multiply(A[len][len - 1]);
                    if (x.getSumAbs() > 0.0) {
                        final Complex y = (A[len - 1][len - 1].subtract(S)).divide(2);
                        Complex z = ((y.multiply(y)).add(x)).getSqrtPositive();
                        if (((y.getRealPart() * z.getRealPart()) + (y.getImaginaryPart() * z.getImaginaryPart())) < 0.0)
                            z = z.getNegation();
                        S = S.subtract(x.divide(y.add(z)));
                    } // end if;
                } // end if; // 400
                for (int i = 0; i <= len; i++)
                    A[i][i] = A[i][i].subtract(S);
                T = T.add(S);
                ++its;
                --sizeMul30;

                // look for two consecutive small sub-diagonal elements
                double xr = A[len - 1][len - 1].getSumAbs();
                double yr = A[len][len - 1].getSumAbs();
                double zr = A[len][len].getSumAbs();
                int m = k;
                // for mm in reverse j..ien-1 loop // 460
                for (int mm = len - 1; mm >= (k + 1); mm--) {
                    print("460 mm=" + mm + ", m=" + m + ", k=" + k);
                    final double yi = yr;
                    yr = A[mm][mm - 1].getSumAbs();
                    final double xi = zr;
                    zr = xr;
                    xr = A[mm - 1][mm - 1].getSumAbs();
                    if (yr <= (((ACCURACY * zr) / yi) * (zr + xr + xi))) {
                        m = mm;
                        break;
                    }
                } // end loop; // 460

                // triangular decomposition A = L*R
                // for i in m+1..ien loop // 620
                for (int i = m + 1; i <= len; i++) {
                    print("620 m=" + m + ", i=" + i);
                    final Complex x = A[i - 1][i - 1];
                    final Complex y = A[i][i - 1];
                    Complex z = null;
                    if (x.getSumAbs() >= y.getSumAbs()) {
                        z = y.divide(x);
                        eigenValues[i] = new Complex(-1.0);
                    } else {
                        // interchange rows of A
                        // for j in i-1..n loop // 540
                        for (int jj = i - 1; jj < size; jj++) {
                            z = A[i - 1][jj];
                            A[i - 1][jj] = A[i][jj];
                            A[i][jj] = z;
                        } // end loop; // 540
                        z = x.divide(y);
                        eigenValues[i] = Complex.get10();
                    } // end if;
                    A[i][i - 1] = z;
                    for (int jj = i; jj < size; jj++)
                        A[i][jj] = A[i][jj].subtract(z.multiply(A[i - 1][jj]));
                } // end loop; // 620

                // composition R*L = H
                // for j in m+1..ien loop // 760
                for (int jj = m + 1; jj <= len; jj++) {
                    final Complex x = A[jj][jj - 1];
                    A[jj][jj - 1] = Complex.get00();

                    // interchange columns of A and vec if necessary
                    if (eigenValues[jj].getRealPart() > 0.0) {
                        // for i in low .. j loop // 660
                        for (int i = 0; i <= jj; i++)
                            ArrayTools.swapArrayIndex1Index2(A[i], jj, jj - 1);
                        // for i in low .. N loop // 680
                        for (int i = 0; i < size; i++)
                            ArrayTools.swapArrayIndex1Index2(eigenVectors[i], jj, jj - 1);
                    } // end if

                    // end interchange columns
                    for (int i = 0; i <= jj; i++)
                        A[i][jj - 1] = A[i][jj - 1].add(x.multiply(A[i][jj]));
                    for (int i = 0; i < size; i++)
                        eigenVectors[i][jj - 1] = eigenVectors[i][jj - 1].add(x.multiply(eigenVectors[i][jj]));

                    // end accumulate transformations
                } // 760
            } // 280

            // a root found
            eigenValues[len] = A[len][len].add(T);
            --len;
        } // end loop; // 260 while

        // all roots found
        // for i in A'RANGE loop
        for (int i = 0; i < size; i++) {
            anorm = anorm + eigenValues[i].getSumAbs();
            for (int jj = i + 1; jj < size; jj++)
                anorm = anorm + A[i][jj].getSumAbs();
        }
        final Complex accnorm = new Complex(anorm * ACCURACY);
        if ((anorm == 0.0) || (size < 2))
            return new MatrixSolverEigen(eigenValues, eigenVectors); // done

        // back substitute to set up vec of upper triangular form
        // for ien in reverse low+1..N loop
        for (len = size - 1; len > 0; len--) {
            final Complex x = eigenValues[len];
            // for i in reverse low .. ien - 1 loop
            for (int i = len - 1; i >= 0; i--) {
                Complex z = A[i][len];
                for (int jj = i + 1; jj < len; jj++)
                    z = z.add(A[i][jj].multiply(A[jj][len]));
                Complex y = x.subtract(eigenValues[i]);
                if (y.getSumAbs() == 0.0) // fp compare to 0 is ok to prevent divide by 0
                    y = accnorm;
                A[i][len] = z.divide(y);
            }
        }

        // multiply by transformation Matrix to give vec of original full Matrix
        for (int jj = size - 1; jj >= 0; jj--)
            // for i in A'RANGE loop
            for (int i = 0; i < size; i++) {
                Complex z = eigenVectors[i][jj];
                for (int k = 0; k < jj; k++)
                    z = z.add(eigenVectors[i][k].multiply(A[k][jj]));
                eigenVectors[i][jj] = z;
            }

        return new MatrixSolverEigen(eigenValues, eigenVectors);
    } // end cxeig2c

    private static MatrixSolverEigen getEigenValuesFrom2x2MatrixComplex(final MatrixComplex matrixComplex) {
        final List<Complex> eigenValues = new ArrayList<>();
        /*
         * To solve: det(A - y*I) = 0
         * <=> (a00 - y)*(a11 - y) - a01*a10
         * <=> a00*a11 - a11*y - a00*y + y^2 - a01*a10
         * <=> y^2 - (a00 + a11)*y + a00*a11 - a01*a10
         * ( := x^2 + b*x + c = 0 => x1/2 = (-b +- sqrt(b^2 - 4*c))/2
         * ( := x^2 - b*x + c = 0 => x1/2 = (b +- sqrt(b^2 - 4*c))/2
         */
        final Complex a00 = matrixComplex.get(0, 0);
        final Complex a01 = matrixComplex.get(0, 1);
        final Complex a10 = matrixComplex.get(1, 0);
        final Complex a11 = matrixComplex.get(1, 1);
        final Complex b = a00.add(a11); // (a00 + a11)
        final Complex c = a00.multiply(a11).subtract(a01.multiply(a10)); // a00*a11 - a01*a10
        final Complex sqrt = (b.multiply(b)).subtract(c.multiply(4.0)).getSqrtPositive();
        eigenValues.add((b.add(sqrt)).divide(2.0));
        eigenValues.add((b.subtract(sqrt)).divide(2.0));
        // return eigenValues;
        return new MatrixSolverEigen(eigenValues, null); // TODO null is not good? -tbach
    }

    private static void print(final String string) {
        if (IS_PRINT_ACTIVATED)
            System.out.println(string);
    }
}
